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The rotation of a macroscopic quantum fluid is a source 
of fascinating problems. By contrast with a classical 
fluid, for which the equilibrium velocity field corresponds 
to rigid body rotation, a quantum fluid described by 
a macroscopic wave function rotates through the nucle- 
ation of quantized vortices 0, • A vortex is a singular 
point (in 2 dimensions) or line (in 3 dimensions) where 
the density vanishes. Along a contour encircling a vor- 
tex, the circulation of the velocity is quantized in units 
of h/m, where m is the mass of a particle of the fluid. 

Vortices are universal features which appear in many 
macroscopic quantum systems, such as superconductors 
or superfluid liquid helium. Recently, detailed investi- 
gations have been performed on rotating atomic gaseous 
Bose-Einstein condensates. These condensates are usu- 
ally confined in a harmonic potential, with cylindrical 
symmetry around the rotation axis z. Two limiting 
regimes occur depending on the ratio of the rotation 
frequency Q and the trap frequency uo in the xy plane. 
When Q is notably smaller than cj, only one or a few 
vortices are present at equilibrium 0, 3 . When Q ap- 
proaches u; , since the centrifugal force nearly balances 
the trapping force, the radius of the rotating gas increases 
and tends to infinity, and the number of vortices in the 
condensate diverges |E El 13 • 

As pointed out by several authors, the fast rotation 
regime presents a strong analogy with Quantum Hall 
physics. Indeed the one-body hamiltonian written in 
the rotating frame is similar to that of a charged par- 
ticle in a uniform magnetic field. Therefore the ground 
energy level is macroscopically degenerate, as the cele- 
brated Landau levels obtained for the quantum motion of 
a charge in a magnetic field. There are two aspects in this 
connection with Quantum Hall Physics. Firstly, when 
the number of vortices inside the fluid remains small com- 
pared to the number N of atoms, we expect that the 
ground state of the system will correspond to a Bose- 
Einstein condensate, described by a macroscopic wave 
function tp(r). This situation has been referred to as 
'mean field Quantum Hall regime' 0, GH [TJ E GH El 



Secondly, when Q approaches uo even closer, the num- 
ber of vortices reaches values comparable to the total 
number of atoms N . The description by a single macro- 
scopic wave function then breaks down, and one ex- 
pects a strongly correlated ground state, such as that 
of an electron gas in the fractional quantum Hall regime 
pjj [3 LL3, [lE Oil • We do not address the second situa- 
tion in this paper and we rather focus on the first regime. 
Furthermore we restrict our analysis to the case of a two- 
dimensional gas in the xy plane, assuming that a strong 
confinement along the z direction so that the correspond- 
ing degree of freedom is frozen. 

The main features of the vortex assembly equilibrium 
in the fast rotation regime are well known. The vortices 
form a triangular Abrikosov lattice in the xy plane and 
the area of the elementary cell is A = irh/(mQ) poj . 
The atomic velocity field obtained by a coarse-grained 
average over a few elementary cells is equal to the rigid 
body rotation result v = H x r, where ft = Qz (z is the 
unit vector along the z axis). 

Beyond this approximation, the physics is very rich 
and many points are still debated. In a seminal paper 0, 
Ho introduced the description of the macroscopic state 
of the rotating gas in the xy plane by a wave function 
belonging to the Lowest Landau Level (LLL). An LLL 
wave function is entirely determined (up to a global phase 
factor) by the location of vortices. Ho considered the 
case of a uniform infinite vortex lattice and inferred that 
the ground state of the system corresponds to a gaussian 
shape for the coarse-grained atom density profile. Using 
both analytical [Til . fl2L I~l| and numerical EH investi- 
gations, it was subsequently pointed out that the atom 
distribution may have the shape of an inverted parabola, 
instead of Ho's gaussian result. Two paths have been 
proposed to explain the emergence of such non gaussian 
profiles. The first one assumes that the restriction to the 
LLL is not sufficient and the contamination of the ground 
state wave function by other Landau levels is responsible 
for the transition from a gaussian to an inverted parabola 
0,0. The second path explores the influence of distor- 
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FIG. 1: Structure of the ground state of a rotating Bose- 
Einstein condensate described by a LLL wavefunction (A = 
3000). (a) Vortex location; (b) Atomic density profile (with 
a larger scale); The reduced energy defined in Eq. (|11.6|> is 
e = 31.4107101. 



tions of the vortex lattice (within the LLL) to account for 
the deviation of the equilibrium profile from a gaussian 

mna. 

In the present paper, we investigate the structure of the 
vortex lattice for a fast rotating condensate. We derive 
the condition under which the LLL is a proper variational 
space to determine the ground state wave function within 
a good approximation. We present a numerical and an- 
alytical analysis of the structure of the vortex lattice, 
based on a minimization of the Gross-Pitaevskii energy 
functional within the LLL. We find that the vortices lie 
in a bounded domain, and that the lattice is strongly 
distorted on the edges of the domain. This leads to a 
breakdown of the rigid body rotation hypothesis which, 
as said above, would correspond to a uniform infinite lat- 
tice with a prescribed volume of the cell. The distortion 
of the vortex lattice is such that, in a harmonic potential, 
the coarse-grained average of the atomic density varies as 
an inverted parabola over the region where it takes sig- 
nificant values (Thomas- Fermi distribution). A similar 
conclusion has also been reached recently in [13, Q . In 
addition to the atomic density profile, our numerical com- 
putations give access to the exact location of the zeroes 
of the wave function, i.e. the vortices. 

An example of relevant vortex and atom distributions 
is shown in Fig. QJi and Fig. ^ for n = 52 vortices. The 
parameters used to obtain this vortex structure corre- 
spond to a quasi-two dimensional gas of 1000 rubidium 
atoms, rotating in the xy plane at a frequency Q = 0.99 cj, 
and strongly confined along the z axis with a trapping 
frequency uj z /(2tt) = 150 Hz. The spatial distribution 
of vortices corresponds to the triangular Abrikosov lat- 
tice only around the center of the condensate: there are 
about 30 vortices on the quasi-regular part of the lattice 
and they lie in the region where the atomic density is sig- 
nificant: these are the only ones seen in the density pro- 
file of Fig.^3. At the edge of the condensate, the atomic 
density is reduced with respect to the central density, the 
vortex surface density drops down, and the vortex lattice 



is strongly distorted. Our analytical approach allows to 
justify this distortion and its relationship with the decay 
of the solution. 

The paper is organized as follows. We start (§U with a 
short review of the energy levels of a single, harmonically 
trapped particle in a rotating frame, and we give the ex- 
pression of the Landau levels for the problem of interest. 
Then, we consider the problem of an interacting gas in 
rotation, and we derive the condition for this gas to be 
well described by an LLL wave function (S ITT)). Sections 
IIIII and ITVl contain the main original results of the paper. 
In § IIIII we explain how to improve the determination of 
the ground state energy by relaxing the hypothesis of an 
infinite regular lattice. We present analytical estimates 
for an LLL wave function with a distorted vortex lattice, 
and we show that these estimates are in excellent agree- 
ment with the results of the numerical approach. In § 11 VI 
we extend the method to non harmonic confinement, with 
the example of a quadratic+quartic potential. Finally we 
give in section IVl some conclusions and perspectives. 



I. SINGLE PARTICLE PHYSICS IN A 
ROTATING FRAME 

In this section, we briefly review the main results con- 
cerning the energy levels of a single particle confined in 
a two-dimensional isotropic harmonic potential of fre- 
quency ou in the xy plane. We are interested here in 
the energy level structure in the frame rotating at angu- 
lar frequency fi (> 0) around the z axis, perpendicular 
to the xy plane. 

In the following, we choose cj, huj, and ^JhJ(rruJ)^ as 
units of frequency, energy and length, respectively. The 
hamiltonian of the particle is 

= -I(V-iA) 2 + (l-tt 2 )y (1.1) 

with r 2 = x 2 +y 2 and A = flxr. This energy is the sum 
of three terms: kinetic energy, potential energy r 2 /2, and 
'rotation energy' —QL Z corresponding to the passage in 
the rotating frame. The operator L z = i(yd x — xd y ) is 
the z component of the angular momentum. 



A. The Landau level structure 

Eq. (ft. 1)1 is formally identical to the hamiltonian of a 
particle of charge 1 placed in a uniform magnetic field 
2Qz, and confined in a potential with a spring constant 
1 — Q 2 . A common eigenbasis of L z and H is the set of 
(not normalized) Hermite functions: 

</>j, k (r) = e r2 / 2 (d x + id y y (d x - id y ) k (e- 1 " 2 ) (1.2) 
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FIG. 2: Single particle spectrum for Q = 0.9. The index k 
labels the Landau levels. 

where j and k are non- negative integers. The eigenvalues 
are j — k for L z and 

E jik = 1 + (1 - Q)j + (1 + Q)k (1.3) 

for H. For Q = 1, these energy levels group in series of 
states with a given /c, corresponding to the well known 
Landau levels. Each Landau level has an infinite degener- 
acy. For Q slightly smaller than 1, this structure in terms 
of Landau levels labeled by the index k remains relevant, 
as shown in Fig. The lowest energy states of two ad- 
jacent Landau levels are separated by ~ 2, whereas the 
distance between two adjacent states in a given Landau 
level is 1 - Q <C 1. 

It is clear from these considerations that the rotation 
frequency Q must be chosen smaller than the trapping 
frequency in the xy plane, i.e. uj = 1 with our choice of 
units. Otherwise the single particle spectrum Eq. (|I.3|) 
is not bounded from below. Physically, this corresponds 
to the requirement that the centrifugal force mQ 2 r must 
not exceed the restoring force in the xy plane —muo 2 r. 



In the LLL, there is a one-to-one correspondence be- 
tween atom and vortex distributions. This relation can 
be made explicit by introducing the atom density p a (r) = 
mr)\ 2 : 

ln(p a (r)) = -r 2 +2^1n|r-r i | . (1.6) 

j 

Introducing the vortex density p v (r) = J2j ~ r j) we 
obtain using V 2 [In \ r — vq\] = 2tt S(r — ro): 

V 2 [ln(p a (r))] = -4 + 4^(r) . (1.7) 

This relation was initially derived by Ho in |9| who in- 
terpreted it in terms of the Gauss law for a system of 
two-dimension charges located at the points rj. 



II. THE INTERACTING GAS IN ROTATION 

We now consider a gas of N identical bosonic atoms 
with mass m. The gas is confined in a cylindrically sym- 
metric harmonic potential, with frequency uj in the xy 
plane and uj z along the z direction. We suppose that the 
characteristic energy huj z is very large compared to all 
other energy scales appearing in the paper, so that we can 
assume that the atoms occupy the ground stat e of the z 
motion, of energy Huj z /2 and extension a z = y/hj (muj z ). 
We are interested in the ground state of this quasi-two 
dimensional gas, when it is rotating at frequency ft close 
to uj around the z axis. 



A. The Gross-Pitaevskii energy functional 



B. The lowest Landau level 

When the rotation frequency ft is close to 1, the states 
of interest at low temperature are essentially those asso- 
ciated with k = 0, e.g. the lowest Landau level (LLL) 
E3 EH E3 Any function i/j(r) of the LLL is a linear 
combination of the </>j,o's and it can be cast in the form: 



^(r) = e- r2/2 P(u) 



(1.4) 



where r — (x,y), u = x + iy and P(u) is a polynomial (or 
an analytic function) of u. When P(u) is a polynomial 
of degree n, an alternative form of ^(r) is 



ip(r) 



2 /2 



Y[(u-Uj) 



(1.5) 



where the Uj (j = 1 . . . n) are the n complex zeroes of 
P(u). Each Uj is the position of a single-charged, positive 
vortex, since the phase of ^(r) changes by 2tt along a 
closed contour encircling Uj. 



The state of the gas is described by a macroscopic wave 
function ip(r) normalized to unity, which minimizes the 
Gross-Pitaevskii energy functional. We introduce the di- 
mensionless coefficient G characterizing the strength of 
atomic interactions, proportional to the atom scattering 
length a s : G = \ /r 87rNa s /a z . The average energy per 
atom, written in the frame rotating at frequency ^, is: 



m = [r 



where is defined in Eq. (|I.1|) . The wave function 

i/)(r) minimizing E[ip] satisfies the Gross-Pitaevskii equa- 
tion: 



d 2 r 



(II-l) 



H^^(r) + G\i>(r)\ 2 ^(r) = ^i,(r) 



(11.2) 



The chemical potential \i is determined by imposing that 
/ = 1. The solution of Eq. (ITL2)) depend s on the two 
independent dimensionless parameters Q and G. 
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The LLL limit 



and tp(r) = outside. The reduced energy is 



In the presence of repulsive interactions (G > 0), the 
basis of Eq. (ft. 2)1 is not an eigenbasis of the TV-body 
hamiltonian. However for a given interaction strength 
G and for a sufficiently fast rotation (Q close to w = 1), 
the restriction to the LLL is sufficient to determine with 
a good accuracy the ground state of the system and its 
energy. Indeed when Q approaches uj = 1, the centrifu- 
gal force mQ 2 r nearly compensates the trapping force 
—muo 2 r and the area occupied by the atoms increases. 
The effect of interactions gets smaller so that the total 
energy per particle tends to the energy Hlj of the lowest 
Landau level, i.e. 1 in our reduced units. 

When ip is chosen in the LLL, the energy functional 
Eq. ()IL1)I can be notably simplified. Indeed the LLL 
functions satisfy the equalities: 



(^kin) — (Eho) 



^ [L z ip] d 2 r 



(11.3) 



where the kinetic and harmonic oscillator energies are: 

<£ki„> = \ J |W| 2 d 2 r (E ho ) = \jr 2 H 2 d 2 r . 

(II.4) 

The total energy E[ip] = is then given by 

E LLL m-Q = J ((l-ft)r 2 H 2 + f M 4 ) d 2 r. (11.5) 

In section IIIII we will minimize this energy functional 
for functions in the LLL. Here we simply outline some 
relevant scaling laws in this regime. 

C. Scaling laws and lower bound in the LLL 

The minimization of Eq. (ftl.5|) is equivalent to the min- 
imization of the reduced energy 



with 



r 2 H 2 + ^H 4 ) d 2 r (II. 6) 



A: 



G 



1-0 



(II.7) 



Therefore the minimizer ^lll depends only on the pa- 
rameter A. This is quite different from what happens 
when the LLL limit is not reached: for the minimization 
of Eq. (ftl.l)) . the two parameters G and Vt are relevant, 
and not only their combination A. 

It is instructive to consider the minimum of e[ip] when 
i/j is allowed to explore the whole function space of nor- 
malized functions J \ip\ 2 d 2 r = 1. This minimum is 
straightforwardly obtained for \ip\ 2 varying as an inverted 
parabola in the disk of radius Rq: 



l^min(r)| S 



ttR 2 



-52 ) ' ^0 - ( — 



1/4 



(II. E 



2\/2 



(11.9) 



The variation of the atomic density as an inverted 
parabola is very reminiscent of the Thomas- Fermi distri- 
bution for a condensate at rest in a harmonic potential. 
However this analogy should be taken with care. In the 
usual Thomas-Fermi approach, one neglects the kinetic 
energy term and the equilibrium distribution is found as 
a balance between potential and interaction energies. In 
the LLL problem considered here, kinetic and potential 
energies are equal (see Eq. ((II. 3)) ). and their sum j r 2 |^| 2 , 
which is large compared to 1 when A ^> 1, is nearly bal- 
anced by the rotation term — Q(L Z ). 

The function ^ m i n clearly does not belong to the LLL, 
since the only LLL function depending solely on the ra- 
dial variable is exp(— r 2 /2). Consequently the reduced 
energy Eq. (|II.9|) is strictly lower than the result of the 
minimization of e[ip] with ip varying only in the LLL. 
In other words, the minimization of Eq. ()IL5)) that we 
perform in the next section, amounts to find the LLL 
function which is "the most similar" to ^min, so that its 
reduced energy is the closest to e m i n . For A ^> 1, we shall 
see that e ~ ay/A, where a is a coefficient of order unity 
to be determined. 



D. Validity of the LLL approximation 

Since e ~ aV^A, the ground state energy E of the fast 
rotating gas determ ined within the LLL approximation 
is Q + a<\/ G(l — Q). Therefore the restriction to the LLL 
is valid as long as the excess energy a^J G(l — Q) is small 
compared to the splitting 2hw = 2 between the LLL and 
the first excited Landau level: 

Restriction to LLL if: G(l - ft) < 1 (11.10) 

When this condition is satisfied, the projection of tp on 
the excited Landau levels is negligibly small. 

It is interesting to compare the scaling laws derived 
in the LLL with the exact relations obtained using the 
virial theorem. For a 2D gas, this theorem gives for the 
ground state of the (possibly rotating) system 

(E ho ) = (E kin ) + (E int ) (11.11) 
(Eint) = § J H 4 d 2 r (11.12) 

while we expect for LLL wavefunctions for A » 1 
(E ho ) = (E kin ) ~ y/A > (E int ) ~ y/G{\ - fi) . (11.13) 

Therefore, within the LLL validity domain of Eq. (ftl.10)) . 
the scaling laws for the predicted LLL energies agree with 
the constraints imposed by the virial theorem. 
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III. THE LLL EQUILIBRIUM DISTRIBUTION 

This section is devoted to the minimization of the re- 
duced energy given in Eq. (|II.6|) for LLL wave functions. 
We start with wave functions corresponding to an infinite 
regular vortex lattice and we derive the corresponding 
energy. Then, we give numerical results which we use in 
the rest of the section as a guide to improve our choice 
for trial wave functions and analyze the distortion of the 
lattice. 



A. The case of a regular vortex lattice 

1. The average density profile for a regular vortex lattice 

We consider a wave function in the LLL with an infinite 
number of vortices on a regular lattice and an average 
spatial density p v . We denote by Uj the points of the 
regular triangular lattice, and by A = l/p v the area of 
its unit cell. We consider the LLL wave functions 



UjeD R 



Uj) , 



(III.l) 



where only the i^ 's located in the disk Dr of radius R 
centered at the origin contribute to the product and the 
constant C is due to the normalization J \yj>\ 2 = 1. For 
A > 7r, we now prove the following result for the atomic 
density p a (r) = |^(r)| 2 : 



Pair) -> p(r) pair) as R 



oo , 



(III.2) 



where p(r) is periodic over the lattice and vanishes at the 
lij's, and where 



Pair) 



1 



-r 2 /* 2 



7T 

A 



(III.3) 



The function p a {r) is the coarse-grained average of the 
atomic density pair). This gaussian decay has already 
been obtained by Ho in the so-called averaged vortex ap- 
proximation Q. However we find useful to prove it here 
with a different approach, which we shall generalize to 
non uniform lattices (§ IIII C)L 

To prove this result, we write ln(p a (r)) = v(r) +w(r), 
with w(r) = ln(p a (r)) and 



\u — uA — ^- [ In \u — u'\ d 2 r r 
AJ Pr 



v(r) = i + 2 22 ln| 

ujeD R 

voir) = -7-r 2 + 4 / ^\u-u f \d 2 r f , (III.4) 



A 



Pr 



where we set u' = x' + iy', 7 = m(7r<r 2 ) and 7' = 7 + lnC. 
Here Pr denotes the inner surface of the polygon formed 
by the union of all elementary cells having their center 
uj in the disk Dr. We want to find the limit of v(r) and 
w(r) as R is large. 



We start with the calculation of v(r). The integral 
entering in the definition of v can be written 



1 

A 



/ ln\u-u f \ d 2 r' 



Pr 



E 

ujeD R 



In \u 



(III.5) 



where the sign f stands for the integration over the unit 
cell of the lattice divided by the area of the cell A. When 
R tends to infinity, v(r) tends to the series 



Vooir) = 70 + 



\U — Uj 



d 2 r v 



U — Uj 



(III.6) 



whose convergence can be checked by expanding the func- 



tion In lu- 



mp to third order in u" / (u - 



This 



series is a periodic function over the lattice and we set 
p(r) = exp(v 00 (r)), which is also periodic. 

To calculate w(r), we first consider the auxiliary func- 
tion w(r) = voir) — wi0)+r 2 /a 2 . Using V 2 [In \r — vq\] = 
27T(5(r — ro), we find that w is harmonic in Pr, with 
w(0) = 0. Moreover, a small computation leads to the 
inequality w(r) > — 7rr 2 /(2*4). In the limit R — > 00, 
we find that w converges to Woo, which is a harmonic 
polynomial with degree less than 2. Due to the sym- 
metry properties of the unit cell, and the lower bound 
by the parabola — 7rr 2 /(2.4), = 0, hence the result 
Eq. (ITTO) . 

To summarize, when the vortex lattice is periodic with 
a uniform average spatial density p Vl the coarse-grain 
average p a of the atomic density is the gaussian of width 
a. The relation in (Eq. 1111.3)1 can be put in the form 



V 2 [ln(p a (r))] = -4 + 47^ 



(III.7) 



which generalizes to coarse-grained quantities the result 
given in Eq. ()L7)I . The fast rotation limit corresponds to 
the case of a large spatial extent of the atom distribution, 
i.e. <j — > +00 or equivalently A = l/p v —> 



2. The energy associated with a uniform vortex lattice 

Once the behavior of the limiting function yj is known, 
we can determine the reduced energy (|II.6)) in the limit 
of fast rotation. This requires the calculation of the in- 
tegrals J pair) d 2 r, j r 2 pair) d 2 r and J p 2 (r) d 2 r in 



the limit R —> 00. It is performed using Eq. (|III.2f) and 
(|III.3)h by taking advantage on the difference in the scales 
of variations of pair) (scale a ^> 1) and p(r) (scale ~ 1). 
We get 



Pair) d 2 r - (fp(r) d 2 r^j x (^J p a {r) d 2 r^j (III. 



8) 



so that the normalization of p a entails f p(r) d 2 r = 1. A 
similar splitting between p and p a occurs for the energy 
and we find: 



J (r 2 Pair) + yPaW) d 2 r = a 



6 



where we have set 



p 2 (r) d 2 r 



(111.10) 



The reduced energy Eq. (IiTL9|) depends on the area A 
of the unit cell through a and on its shape through the 
Abrikosov coefficient b. Let us briefly recall the origin of 
this coefficient. Instead of using the exact atomic density 
p a (r) to calculate the energy, we work with the coarse- 
grain average p a (r), whose spatial variation is much sim- 
pler. To do this substitution, we must renormalize the 
interaction coefficient G, which becomes bG. This is a 
consequence of the discreteness of the vortex distribu- 
tion: since the wave function ip(r) must vanish at the 
vortex location, the average value of |^| 4 over the unit 
cell, hence the interaction energy, is larger than the result 
obtained if was quasi-uniform over the cell. 

We now look for the choice of b and a which minimizes 
the reduced energy Eq. QUI. 9)1 . As known for the case 
of superconductors, the lattice minimizing b is the trian- 
gular one 24], for which b ~ 1.1596. The minimization 
over a then leads to: 

do = (6A/(4tt)) 1/4 and e = y/bA/ir. (III.ll) 

We recover a scaling similar to Eqs. (|II.8)) - f|II.9)) . inferred 
for a distribution varying as an inverted parabola. Note 
that the size of the elementary cell A = 7r(l — <Tq~ 2 ) -1 
differs from the rigid body rotation result, *4rbr = 7r/f2, 
although the two quantities tend to tt when Q tends to 1. 
Actually if we impose A = *4rbr in Eq. QUI. 9)) . instead 
of minimizing on <j, we find that 1/cr 2 = 1 — f2 and we 
obtain ~ 2, much larger than the result ~ 1 

deduced from Eq. QUI. 11)1 . 

The reduced energy eo exceeds the lower bound 
Eq. (ITL9)) by the factor Vb x ^/9/S - 1.14. The ori- 
gin of the c oeffic ient b has been explained above. The 
coefficient y/9/8 = 1.06 is due to the difference between 
the gaussian envelope found here (c.f. Eq. ((111.3)1 ). and 
the optimum function varying as an inverted parabola 
Eq. (FTTH)) . For the parame ter A 3000 used in Fig. [TJ 
we find eo = 33.3 using Eq. QUI. 11)1 . which is ~ 6% larger 
than the result found numerically (cf. Fig. ^ an d § IIII Bl 
below). 



B. Minimization in the LLL: numerical results 

We now turn to the description of the numerical 
method that has been used to obtain the vortex and 
atomic patterns shown in Fig. E] and we give some fur- 
ther results of interest for the following discussion. For a 
given A = G/(l — Q) and a given number n of vortices, 
we write our trial functions under the form of Eq. ((1.5)) . 
We vary the location of the vortices U{ using a conjugate 
gradient method to determine the optimal location and 
the minimum reduced energy e n . 

The computation of the energy uses the Gauss point 
method (see 25] or 26] for a use in the case of the 
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FIG. 3: Minimum reduced energy e n as a function of the 
number of vortices in the trial wave function (A = 3000). 



Gross Pitaevski equation): the computation of the in- 
tegral of a polynomial times a gaussian is exact as long 
as the degree of the polynomial is lower than a certain 
bound, which depends on the number of Gauss points. 
An alternative method used for example in [l4( consists 
in writing the trial functions in the form Eq. ((1.4)1 with 
P( u ) — Y^j=i bjui , and performing the minimization by 
varying the coefficients bj . The advantage of the method 
followed here is to give directly the location of the vor- 
tices, while the alternative approach requires to find the 
n roots of the polynomial P(u), which may be a delicate 
task for large n. 

For the range of A's that we have explored (between 
300 and 3000), the reduced energy e n decreases for in- 
creasing n, until it reaches a plateau. For A = 3000 
(Fig. |3J), the plateau is reached for n = 52 and the re- 
duced energy varies in relative value by ~ ±10 -8 when 
n increases from 52 to 70. The vortex and atom distri- 
butions for n = 52 are given in Fig. Q When n increases 
the central distribution of vortices remains the same, as 
well as the significant part of the atom distribution. The 
distribution minimizing the energy for n = 70 vortices 
is shown in Fig. We note that beyond n = 52, the 
location of the additional vortices strongly depends on 
the initial data of the optimization procedure, as extra 
vortices only change slightly the energy. In addition to 
the result of Fig. which is the absolute minimum for 
{A = 3000, n = 70}, we have found a number of config- 
urations corresponding to local minima where the addi- 
tional vortices lie on an outer distorted circle. 



C. The distorted lattice 

Inspired by the numerical results such as the ones 
shown in Fig. and Fig. QJ we generalize the approach 
developed for the regular vortex lattice to the case of 
a distorted lattice. We make the hypothesis that the 
locations Uj of the vortices are deduced from a regular 
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FIG. 4: Vortex distribution minimizing the reduced energy 
for A = 3000 and n — 70 vortices. 

hexagonal lattice u™ s by 

Uj = r(«f g ) = A(|«f e |)«f g (ni.12) 

where A(r) is a positive function varying smoothly over 
a distance of order unity. We assume that the unit cell 
of the initial regular lattice has the area A = 7r, corre- 
sponding to a flat density profile in Eq. QUI. 3)1 . If A tends 
to infinity for a finite value r^, the number of vortices in 
the distorted lattice is finite and equal to ~ r^, since all 
vortices located after the 'horizon' rh in the regular lat- 
tice are rejected to infinity. Otherwise, if A is finite for 
all r, the number of vortices in the distorted lattice is 
infinite. 

The distortion is illustrated in Fig. [5] for the particular 
case of A = 3000. We have plotted at the same scale 
the regular lattice with A = tt and the configuration 
of vortices minimizing the energy. For n = 52 vortices, 
only the lattice sites of the regular lattice whose distance 
to the origin is below rh = 7.4 remain in the distorted 
lattice. Around the center of the disk of radius r^, the 
function A(r) is close to 1, whereas it becomes very large 
when r approaches the 'horizon' rv 

As for the case of the regular lattice, we introduce the 
coarse-grained averages p a and p v of the atom and vortex 
densities. The function p v is now space dependent and 
is simply the inverse of the area of a distorted cell in the 
vicinity of r : 

p v (r) = (7rA(r')(A(r') + r'A'(r')))" 1 , (m.13) 

where r = \{r')r f . We recall that the expected length 
scale in the limit of fast rotation is Rq = (2A/7T) 1 / 4 ^> 1 
(see § III C)l and we consider a class of distortions A(r) 
such that 

A 2 (r)=l + ^^ + o(^) , (111.14) 

where f((, 2 ) is a continuous function, which diverges at 
£h = r h/Ro- We also assume that the integral F(s) = 
Jo f( s ') ds' diverges at s = We shall check in the end 
that the distortion minimizing the energy belongs to the 
class of functions defined in Eq. QUI. 14)1 . 
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FIG. 5: Regular lattice with A = tv and distorted lattice 
minimizing the energy for A = 3000 and n — 52 vortices. 

In the limit A ^> 1, we show in the appendix the fol- 
lowing properties for the vortex lattices obtained through 
a distortion obeying (Eq. 1111.14)) : 

1. The atom density p a (r) can be written as 

hi(p a (r))=v(r)+w(r) (111.15) 

where v is related to the function Voo(r) introduced 
for a regular lattice in Eq. QUI. 6)) : 

v(r) = Voo(r') with r = \(r')r' . (111.16) 

w(r) is a smooth radial function and we set p a (r) = 
exp(w(r)), where p a is normalized to unity. 

2. The coarse-grain average p a at a point r = Ro£ 
and the integral F of the distortion function / are 
related by the relation: 

p a (RoO*eM-F(t 2 )) if (IIL17) 

and is zero elsewhere. Note that p a (r) is continuous 
at rh = RoCh since we have assumed that F(s) 
tends to +oo when s — > 

3. As for the regular lattice case, we use the difference 
in the scales of variations of the two functions v and 
w to obtain 

J (^ 2 p a (r) + b -±pl(r))d 2 r. (111.18) 

where b = 1.1596 as for a regular lattice. 

The differences with respect to the initial minimiza- 
tion problem of Eq. (|II.6j) are the renormalization of the 
coefficient G — > bG discussed in § IIII Al and the fact that 
p a is a smooth, non-negative radial function, instead of 
being the square of an LLL wavefunction. 

D. The "Thomas-Fermi" distribution in the LLL 

We now address the minimization of the energy func- 
tional in Eq. 1)111. 180 . The minimizing function is the in- 
verted parabola p a (r) ex R\-r 2 for r < R 1 = (26A/7T) 1 / 4 , 
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and p a = for r > R\. The associated energy is 

e c± — = 
3v^ 



6A. 



(111.19) 



Using Eq. fflll. 1 T|) we deduce the distortion function f(s) 
and its primitive F(s): 



f(s) 



Vb- 



F(s) 



In 1 



Vb 



(111.20) 



As initially assumed, the functions f(s) and F(s) tend 
to +oo at the horizon y/b, hence rh = b^^R^ = R\. This 
means that at leading order in A, the Thomas Fermi 
radius and the horizon are equal. The function T trans- 
forming the initial regular lattice u^ g into the distorted 



lattice Uj is thus: 




FIG. 6: Radial density distribution p T &d(r) for A = 3000. 
The unit along the vertical direction is arbitrary. The dotted 
line is a fit using the inverted parabola with the radius Ri = 
(26A/7r) 1//4 , with an adjustable amplitude. 



r = T(r') = r' 



R(-r 



,/2 * 



(111.21) 



Once f(s) is known, one can evaluate the vortex density 
using Eq. (1111. 130 : 



Ri 



(Ri 



~/2\2 



(111.22) 



where r and r' are related by Eq. (|III.21|) . In particular 
the vortex density at r = is ~ (1 — R~^ 2 )/ii, which is 
close, but not equal, to the prediction Q/tt for a rigid 
body rotation. 

Our distortion function / is to be related to that of [2^ , 
though it is derived using very different techniques. The 
asymptotic result Eq. (|III.19|) has also been obtained re- 
cently by Watanabe, Baym and Pethick 13] who assumed 
that Eq. (|TlL7|) can be generalized to the case where p v 
is spatially dependent: 



V 2 [ln(p a (r))] =-4 + 4^(r) 



(111.23) 



By differentiating Eq. (|III.17|) . a similar relation can 
be proved within our approach with p v (r) replaced by 
p v (T(r)). The two relations are equivalent at points not 
too close to the Thomas- Fermi radius (i.e. R\ — r > 1). 
A result related to Eq. (|III.23|) has also been shown in a 
different context by Sheehy and Radzihovsky |27j |. They 
consider the case of a condensate which is not in very fast 
rotation (i.e. outside of the LLL regime) but still with 
several vortices. Interestingly, the procedure used in [2^ 
to derive the relation between p v and p a is based on the 
minimization of the energy functional, including atom 
interactions. On the contrary, the result in Eq. (|III.17|) 
or Eq. (|III.23|) is a consequence of the structure of an 
LLL wave function and it is at first sight independent 
of atomic interactions. However one must keep in mind 
that the knowledge of the strength of atom interactions 
is essential to check the relevance of LLL wave functions 
for the problem (see Eq. (11.10)). The relation reached 
in [2jj has the same structure as Eq. (|III.23|) . but with a 



dimensionless coefficient involving the healing length and 
p v inside the V 2 ln(p a ) term. Close to the Thomas- Fermi 
radius, p v varies rapidly and the approach of [2jj leads to 
a different relation from Eq. (|III.23|) , since the derivatives 
of p v have a significant contribution in this region. 

Our analytical predictions can be compared with our 
numerical results obtained in the particular case A = 
3000 (i.e. R x = 6.86) , for w hich we plotted Fig. □ 
The prediction of Eq. (1111. 19|) yields e = 31.374, only 
0.12% below the value determined numerically. We can 
also compare our trial density with the numerical result. 
We give in Fig. the prediction of the inverted parabola 
together with the radial density distribution determined 
numerically: 



Prad(r) 



2tt 



2tt 



|^(r)| 2 d0 



(111.24) 



where is the polar angle in the xy plane. Apart from os- 
cillations due to the discreteness of vortices, the two dis- 
tributions are remarkably close to each other. A similar 
conclusion was reached recently by Cooper, Komineas, 
and Read [l4|. They also performed a numerical mini- 
mization of the energy of Eq. (|II.1|) in the LLL limit, and 
found an atom density profile in good agreement with 
the inverted parabola distribution predicted in . 

From the above analytical results, we expect that the 
minimizing configuration will involve n ~ r\ ~ 48 vor- 
tices. The number of vortices for which the minimum en- 
ergy plateau is reached numerically is 52, which is very 
close to r 2 . As for the location of vortices, our analy- 
sis indicates that the vortices in the distorted lattice are 
images through Eqs. (|III.12|) - (|III.21|) of points of the reg- 
ular hexagonal lattice such that |i^ eg | < = R\. Note 
that the optimal vortex configuration involves some vor- 
tices outside the disk of radius R\. They correspond to 
regular lattice sites |^ eg | close to the horizon r^. In- 
deed, for these points A gets large and the image point 
is sent beyond the Thomas Fermi radius. Thus, though 
the distorted lattice provides an inverted parabola which 
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vanishes at Ri, the location of the vortices extends be- 
yond R\. The numerical analysis leads to results which 
nicely confirm our analytical predictions. In addition it 
allows to explore the role of the vortices lying outside 
the Thomas-Fermi distribution. For example one can re- 
move the contribution (u — Uj) of these vortices in the 
expression Eq. (ft. 5)) of the LLL wave function, while keep- 
ing unchanged the contribution of the vortices inside the 
Thomas- Fermi radius. This results in a significant mod- 
ification of p a (r) which then vanishes around ^ 7.3, in- 
stead of ~ 6.8. Therefore these outer vortices play an 
important role in the equilibrium shape of the conden- 
sate, even though they cannot be found when one simply 
plots the atomic spatial density. 

A closer look at Fig. indicates that p a is matched 
to zero more smoothly than an inverted parabola. 
An expansion of the energy of the distorted lattice 
(Eq. QUI. 18)1 ) to the next order in A should lead to a 
minimizing function p a with a smoother decay to zero 
around R\. In particular a natural way to match the in- 



verted parabola with the asymptotic decay r 



2n, 



of 



any LLL function with n vortices, could be obtained 
through a Painleve-type equation (as at the border of 
a non-rotating BEC). 

Remark: Comparison with the "centrifugal force ap- 
proximation" '. Under some conditions, it is possible to 
write iJj as the product of a rapidly varying function rj(r) 
and a slowly varying envelope ip(r) [ll|. This is rem- 
iniscent of the splitting of ln(p a ) in terms of v and w, 
although it leads to a different conclusion. One obtains 
for the envelope an equation similar to Eq. (ftl.2)h where 
only the centrifugal potential remains jl l| : 

_I V 2 ^(r) + (1 - n 2 )^(r) + G\4>(r)\ 2 ^(r) = fi$(r) 

(111.25) 

where Jjl = \i — O. We call this approach the "centrifugal 
force" approximation [2^ and we compare its predictions 
with those derived from the LLL approximation. 

When the approximation leading to Eq. QUI. 25)1 is 
valid, one is left with the problem of a 2D gas at rest 
in a harmonic potential with the spring constant 1 — Q 2 . 
The solution of this equation depends on the strength of 
the interaction parameter G. If G ^> 1, the kinetic en- 
ergy term can be neglected (Thomas-Fermi approxima- 
tion) and one gets \i/j(r)\ 2 cx l — r 2 /R 2 fa inside the disk of 

radius i? cfa = (4G/[tt(1 - tt 2 )]) 1/4 and ip(r) = outside. 
Note that i? c f a coincides with our Thomas-Fermi radius 
Ri for Q ~ 1. If G <C 1, the interaction term can be 
neglected and the solution is the ground state of the har- 
monic oscillator, i.e. the gaussian of width (1 — SI 2 ) -1 / 4 . 

In the LLL, we have seen that the distinction between 
the two regimes G ^> 1 and G <C 1 is not relevant. The 
only important parameter is A = G/(l— O). When A ^> 1 
the envelope of the atom density profile is close to an in- 
verted parabola, irrespective of the value of G. Therefore, 
there exists a clear discrepancy between the predictions of 
the LLL treatment and those of the centrifugal force ap- 



proximation when 1 — Q <C G <C 1 . For these parameters 
the LLL approximation is valid since G(l - (]) < 1 (see 
Eq. jnU)- The extent of the wave function minimizing 
e[ijj] is thus Ri ~ (G/(l — SI)) 1 / 4 , while the reasoning 
based on Eq. QUI. 25)) would lead to a gaussian envelope 
with a larger size (1 — ^) -1 / 4 , independent of G. 



IV. EXTENSION TO OTHER CONFINING 
POTENTIALS 

The ideas that we have developed for a harmonic con- 
finement can be generalized to a larger class of trapping 
potential V(r). For simplicity we assume here that V 
is cylindrical symmetric, with a minimum at r = 0. We 
define uo as moo 2 = d 2 V/dr 2 \o and we set 



V(r) = \muo 2 r 2 + W{r) . 



(IV.l) 



As above we choose uo and y/h/(mu) as the units for 
frequency and length, respectively. 

We are still interested here in a region where £1 ~ 1. 
To minimize the Gross-Pitaevskii energy functional, we 
use again wave functions in the LLL so that the energy 
per particle to be minimized is 



LLL 



((l-n)r 2 + W(r))p a + ^p 2 a 



d 2 r 



(IV.2) 

As explained in section II. C, the LLL approximation is 
valid if the minimum for £^lll — ^ is small compared 
to the distance 2 = 2huj between the LLL and the first 
excited Landau level. 

We have seen that varying the locations U{ of the vor- 
tices, hence the average vortex surface density p v , allows 
to generate a large class of coarse-grain averaged atom 
densities p a . Provided W(r) is well behaved, we can gen- 
eralize the treatment presented for the purely quadratic 
case. The energy £^lll can still be expressed in terms of 
p a instead of p a with an expression similar to Eq. (ft V. 2)1 . 
and the interaction parameter G replaced by bG. The 
Thomas-Fermi distribution minimizing ^lll is 



pi F (r) 



max 



/i- (1 - n)r 2 - W(r) 
bG 



(IV.3) 



where p is the chemical potential determined such that 
J p a = 1. Once p a has been determined over the whole 
space, the energy I^lll can then be calculated and the 
validity of the various approximations can be checked: 
(i) I^lll -^| < 1 and (ii) the extension of the domain 
where p a differs from zero is large compared to 1, so that 
it is legitimate to introduce a coarse-grain average of p a 
over several vortex cells, and there is a large parameter 
playing the role of Rq. 

As an example, we investigate the case of a combined 
quartic and harmonic potential: W(r) = /cr 4 /4, which 
has been studied recently both theoretically, numerically 
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HI HI IH |H Hi EI El and experimentally H3- A 
nice feature of this potential is that it allows to explore 
the region ft > 1, since the centrifugal force, — ft 2 r, can 
always be compensated by the trapping force, varying 
as -(r + kr 3 ). We define A = (3/c 2 6G/(8tt)) 2 / 3 and 
A = (1 — ft) 2 + kfjL. Two cases can occur, (i) If ft < 
ft c = 1 + a/Ao? then p^ F is non zero in a disc of radius 
R\ = 2 (ft - 1 + >/A)/A, ft and A being linked by 

2A 3/2 + 3A(ft - 1) - (ft - l) 3 = 4Aq /2 . (IV.4) 

(ii) If ft > ft c , then pJ F is non-zero on an annulus of radii 
R 2 ± = 2(tt-l± Vk)/k, and A = A . 

The Thomas- Fermi distribution given in Eq. (|IV.3|) al- 
lows to calculate the minimum energy per particle. Since 
the general calculation is quite involved, we simply give 
here the result for ft = ft c : 

ft = ft c : E LLL - ft = a k 1/s G 2/3 (IV.5) 

where a ~ —0.1. More generally, when |1 — ft| is at most 
of the order of /c 2/3 G 1/3 , then ^lll - ft is of the order of 
q2/3 rpj^ res t r i c ti on t Q the LLL wave functions and 
the use of the 'Thomas-Fermi' approximation (Eq. (|IV.3|) 
are valid if two conditions are fulfilled: (i) — ft <C 

1, hence kG 2 <C 1, (ii) the extension R + ~ (G/k) 1/6 
of p a is large compared to 1, so that the coarse-grain 
average of p a is meaningful. This requires « 1 and 
k <C G <C 1/Vk. When these conditions are satisfied, 
ft c - 1 - k 2/3 G 1/s <C 1, and the study of the regime 
ft > ft c can be performed within the LLL. In addition 
one can check that for ft c - 1 < ft - 1 <C /c 1/3 G 2/3 , the 
width R + — R_ of the annulus is large compared to 1 
(both R+ and R_ are of order (G/Zc) 1 / 6 ), so that the use 
of the coarse-grain averages of p a and p v is justified. A 
similar analysis to what we have performed above yields 
an almost uniform vortex lattice in the annulus, with a 
distortion near the inner and outer boundaries. 

The LLL approximation has been used by Jackson, 
Kavoulakis and Lundh to study the phase diagram of the 
vortices in a quadratic+quartic phase 35]. They were 
mostly interested in the stability of giant vortices, hence 
they restricted their analysis to particular LLL states, 
where F(u) only contains two or three terms bjuK How- 
ever one could in principle use the same approach as the 
numerical treatment developed here, and derive the de- 
tailed vortex pattern for various choices of G, ft and k. 
It would be interesting to see whether there exists a do- 
main of parameters where the polynomial F(u) has a 
multiple root in u = 0. This would correspond to the gi- 
ant vortex which has been predicted by other approaches 
HO, El H^. Another limit where R + - R_ < 1 has re- 
cently been studied in [36^ . 

V. CONCLUSION 

In this paper, we have studied analytically and nu- 
merically the vortex distribution and atomic density for 



the ground state of a rotating condensate trapped in a 
harmonic potential, when the rotation and trapping fre- 
quencies are close to each other. Restricting our anal- 
ysis to quantum states in the lowest Landau level, we 
have shown that the atomic density varies as an inverted 
parabola over a central region. The vortices form an 
almost regular triangular lattice in this region, but the 
area of the cell differs from the prediction for solid body 
rotation. In the outer region, the lattice is strongly dis- 
torted. We have determined the optimal distortion, and 
related it to the decay of the wave function close to the 
Thomas- Fermi radius. 

Our results agree with those of a recent numerical 
study pjl- Another analytical approach to this prob- 
lem has recently been given in [l3|. It leads to the same 
value as ours for the energy of the ground state, whereas 
our treatment provides more detailed information on the 
vortex pattern at the edge of the condensate. Our pre- 
dictions for the equilibrium shape of the atomic density 
and for the vortex distribution should be experimentally 
testable. In 7] the regime of fast rotation in the LLL 
has already been achieved and it was indeed found that 
the atom density profile varies as an inverted parabola, 
and not as a Gaussian as one would expect for an infinite 
regular lattice |?j • In , a detailed experimental analysis 
of the vortex spacing as a function of the distance to the 
center of the trap has been made and it showed a clear 
distortion of the pattern on the edges of the condensate. 
This study was not performed in conditions such that 
our LLL approximation is valid, and the relevant theo- 
retical model is rather the one developed in J2?| . However 
it should be possible to perform a similar experimental 
analysis for faster rotation rates, and test in particular 
the validity of our prediction concerning the distortion 
factor A(r) (see Eqs. (IIII.12|) - (jlII.14^flllL20^ . 

Finally, we have addressed the case of other trapping 
potentials, such as a superposition of a quadratic and a 
quart ic potential, which have also been addressed experi- 
mentally |3?| . For even faster rotations, when the number 
of vortices approaches the number of atoms, the ground 
state is strongly correlated. We did not touch this point 
here, but our work should be relevant for studying the ap- 
parition of this correlated regime from a destabilization 
of the mean field results by quantum fluctuations. 
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FIG. 7: Example of a distorted lattice generated by the trans- 
formation Eq. (|VI.1J) . The radius of the circle is X a aRo. In 
the regular part outside the circle, the cell area is A a . 



VI. APPENDIX 

The aim of this appendix is to prove the properties 
used in § IIII CI A detailed proof will be given in [38] . We 
consider a distorted lattice in an inner region and keep 
a regular lattice in the outer region in such a way that 
the distortion is continuous (see Fig. [7J). We label by j 
a regular hexagonal lattice with a unit cell area A = tt 
and we define the transformed lattice by 



X(\j\)j for 
X a j for 



< aRo 
> aRo 



(VI.1) 



where the radius Ro is given in Eq. (|II.8|) . the distortion 

and aRo is smaller than 



function A(r) satisfies (|III.14|) 
the horizon where A(r) diverges. We have set X a = 
X(aRo) and the area of the unit cell of the outer lattice is 
A a = ttA^. When aRo tends to r^, A a tends to infinity 
and the vortex lattice of Fig. [7| is similar to the one in 
the right of Fig. [51 

In the following, we shall (i) define p a and compute its 
limit when Ro increases (i.e. Q tends to 1) for a fixed a, 
(ii) let a get close to the horizon £^ = Vh/Ro- We need 
to take the limits in this order, because we will use that 
X a is close to 1, which is only true if a is fixed less than 
£h and Ro is large. 

Firstly we consider only the points j in a disc Dr> 



i/;(r) = Ce~ r ^ 2 H (u 
\3\<R' 



(VI.2) 



Q a denotes the unit cell of the lattice of area A a and 
Pql,R' is the polygon formed by the union of all elemen- 
tary cells of area A a and center X a j, with \j\ < R' . We 



write ln(p (r)) = v R >(r) + w R *(r) with, 

v R ,{r) = 2 (ln|r-A a j| 

aR <\j\<R' 



J_ f ln\r-r'-X a j\d 2 r') 
A* J Qa / 

2 Yl (ln|r-A(j)j| 

\j\<aRo 

j-J \n\r-r' -\{j)j\d 2 r') 



and WRf(r) = w\R>{r) + W2(r) with 



w 2 (r) 



In \r — r'\ d 2 r' 



^ A a Jo \ r - r' - X a i\ 

\j\<ocR a JQa 1 Jl 



We have just added and subtracted terms at this stage. 

Now we let R' tend to infinity and find the limit for an 
infinite number of vortices. This step is very similar to 
the case of the regular lattice since the lattice distortion 
only affects a finite number of sites. We find that 



wiR'(r) - wlr'(O) itfi(r) = -r 2 /a 2 



(VI.3) 



with (j~ 2 = 1 — 7r / 'A a - vr> tends to a convergent series v 
(which is not a periodic function, contrary to the regular 
lattice case). 

The next step is to let Ro be large, keeping a fixed, so 
that X a is close to 1 for the class of distortion functions 
considered in Eq. (|III.14|) . We find 

v(r)^ Voo (r f ) (VIA) 

where Voo is given by Eq. (|III.6|) . and r, r' are related by 



_ J A(r') r' for r/X a < aRo 
X a r' for r/Xa > aRo , 



(VI.5) 



We denote w = + W2 = ln(p a ). We estimate W2(r), 
using an expansion of the logarithm and the fact that 
A(j) ~ X a ~ 1: 

MRot) * - ( (/(« 2 ) - m> 2 )) ^"f . 

TT J S '<a IC-CI 

(VI.6) 

where relevant ^'s are of order unity. Using an integra- 
tion by part and a primitive F of /, we get (6(x) is the 
Heaviside function): 

w 2 {R Q i) ~ [F(a 2 ) - F(C 2 ) + if - a 2 )f(a 2 )] 9{a - 

(VI.7) 

Since we have a 2 ~ f(a 2 )/R 2 ) , then wi(Ro£) = 
—£ > 2 f(a 2 ). Putting everything together, we obtain up 
to an additive constant for normalization, 



-F(f) 



for i < a 



-£ f(a ) + fJL for ^ > a 



(VI.8) 
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with fj, = a 2 f (a 2 )- F(a 2 ). 

Finally, using that X a ~ 1, we can apply the separation 
of integrals (2^ and find for example that 

J p a d 2 r oc (y e ^') dV) x 

The last integral in the second line is equal to 



7re _F( ^ Q;2 ^//(a 2 ). At this stage, a is still a free param- 
eter. If the distortion function A(r) has a horizon at 
r = £fci2o, we let a tend to otherwise to oo. The last 
integral tends to zero given the hypothesis that / and F 
tend to +oo at The same procedure is valid for all 
terms entering into the energy functional, which justifies 
the use of En. fill 18)) . 
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